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I. 


INTRODUCTION 


The  use  of  precision  inertial  sensors  for  high  performance 
applications  demands  accurate  identification  of  drift  coefficients. 

Thic  identification  technique  consists  of  two  parts:  testing  and  data 
reduction.  An  acceptable  identification  technique  should  be  low  cost. 

There  are  at  least  three  different  objectives  for  testing.  First, 
development  tests  are  conducted  in  laboratories  for  verifying  concepts, 
gaining  intuitions,  exploring  desired  modifications,  and  searching  for 
proper  parameter  values.  Second,  acceptance  tests  are  used  by  manufac- 
turers and  customers  to  determine  if  the  specifications  are  met. 

Third,  calibration  tests  are  used  to  determine  sensor  parameters. 

This  report  concerns  test  and  data  reduction  techniques  far 
single-degree-of-f reedom  gyros.  Techniques  for  two-degree-of-freedom 
gyros  are  similar  in  concept  and  therefore  will  not  be  discussed  in 
this  report. 

The  most  important  aspects  of  gyro  testing  for  navigation  and 
guidance  purposes  are  those  concerned  with  drift  errors.  There  are 
at  least  four  different  types  of  testing  methods  for  the  determination 
of  drift  parameters  [1,  2,  3,  4,  5,  and  6].  They  are  as  follows: 

a)  Six  position  rate  test. 

b)  Six  position  static  test. 

c)  Single-axis  tumble  test. 

d)  Two-axis  tumble  test. 


All  methods  require  the  use  of  a test  table.  Both  the  six  position 
rate  and  static  tests  are  time  consuming,  because  they  require  repeated 
reorientation  of  the  gyro  and  a long  settling  time  after  each 
reorientation.  The  single-axis  tumble  test  does  not  generate  sufficient 
information  for  the  determination  of  drift  parameters  of  a high  preci- 
sion drift  model.  Among  the  four  methods,  the  two-axis  tumble  test 
is  the  latest  technique.  Compared  to  the  others,  it  has  several 
advantages  which  will  be  discussed  in  a later  section.  However,  the 
known  two-axis  tumble  test  requires  the  knowledge  of  the  attitude  of 
the  gyro  case  at  every  data  taking  instant. 

This  report  presents  a new  two-axis  tumble  test  and  data  reduction 
technique.  The  method  uses  a new  concept  for  obtaining  all  drift 
parameters  from  the  measurement  record.  The  concept  does  not  need 
the  knowledge  of  the  gyro  case  attitude  at  each  data  taking  instant. 
Thus,  much  simpler  test  equipment  can  be  used.  The  report  will  first 
review  the  development  of  an  analytic  model  for  drift  parameters.  The 
test  procedure  for  the  proposed  new  two-axis  tumble  test  will  be 
described,  followed  by  the  development  of  the  data  processing  algorithm. 
A comparison  between  the  old  and  new  methods  will  be  made.  The  result 
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of  a simulation  using  the  proposed  technique  will  be  included  and  dis- 
cussed. For  the  convenience  of  the  reader,  a list  of  references  and 
an  appendix  are  attached. 


II.  ANALYTIC  MODEL  FOR  GYRO  DRIFT 

Gyro  drifts  may  be  classified  into  three  types:  constant 

drifts,  acceleration  sensitive  drifts,  and  drifts  which  are  sensitive 
to  the  products  of  accelerations  [7]. 

The  constant  drifts  are  caused  by  the  spring  force  of  flex  leads 
and  by  the  electrical  reaction  torques  excited  on  the  gimbal.  All 
constant  drifts  will  be  lumped  together  and  represented  by  Dq. 

Most  of  the  acceleration  sensitive  drifts  are  caused  by  mass 
unbalance  of  the  gimbal-rotor  assembly.  This  unbalance  of  mass  results 
in  a displacement  between  the  gimbal  output  axis  and  the  mass  center 
of  the  gimbal-rotor  assembly.  This  displacement  is  denoted  by  (i  which 
is  a vector  in  the  three-dimensional  space.  Let  s,  i,  and  o be  three 
unit  vectors  pointed  along  the  spin  reference  axis,  input  axis,  and 
output  axis  of  the  gyro,  respectively.  Let  a_  be  the  acceleration 
vector  of  the  gyro  case  and  m be  the  mass  of  the  gimbal-rotor  assembly. 
The  inertial  torque  about  the  output  axis  due  to  the  acceleration  is 
given  by 


T =-m(dxa)*o=m(d.a  -da.) 
m v_  is  Si 


where  "x"  and  denote  cross-product  and  dot-product,  respectively. 
The  subscripts  s and  i are  used  to  denote  those  vector  components  which 
are  in  the  directions  of  s and  i,  respectively.  Because  drift  is 
linearly  proportional  to  torque,  the  mass  unbalance  drift  may  be  repre- 
sented by 


D = K a + K.  a. 
m s s 11 


(1) 


where  Kg  and  are  proportionality  constants. 

Another  source  of  acceleration  sensitive  drift  is  the  thermo- 
convection torque  for  a gyro  with  floated  gimbal.  When  the  mass  is  not 
balanced,  the  thermo-convection  in  the  fluid  generates  a torque  along 
the  output  axis.  This  torque  is  proportional  to  the  mass  m,  off-center 
displacement  d , and  the  acceleration  component  aQ,  where  subscript  o 

A 

denotes  the  vector  components  in  the  direction  of  o.  The  corresponding 
drift  may  be  represented  by 


4 


5 


(4) 


D = D + K a + K.a.  + K a + K a2  + K.  .a2.  + K .a  a . 

o SS  11  OO  SSS  11  1 SI  SI 


+ K.  a. a + K aa 
10  1 o OS  o s 


Equation  (4)  shows  that  the  total  drift  D consists  of  nine  terms, 
characterized  by  drift  parameters  Dq  and  K's.  Because  a different  type 

of  compensation  is  required  for  a different  type  of  drift,  identification 
of  these  parameters  is  essential. 


III.  A NEW  TWO-AXIS  TUMBLE  TEST 

Identification  of  the  nine  drift  parameters  by  testing  requires 
that  all  nine  drift  modes  as  shown  in  Equation  (4)  be  excited  by  test 
signals.  This  can  be  accomplished  by  tumbling  the  gyro  about  two 
orthogonal  axes.  The  proposed  two-axis  tumble  test  consists  of  two 
parts:  performance  of  the  experiment  and  data  processing  for  drift 

parameter  identification.  These  two  parts  are  interrelated,  because 
the  required  measurements  depend  on  the  way  data  are  processed. 

Figure  1 shows  a gyro  test  stand.  The  stand  has  two  rotatable 
axes  which  are  perpendicular  to  each  other.  The  gyro  under  test  is 
mounted  on  the  second  turntable  with  its  input  axis  (IA)  precisely 
aligned  to  the  table  axis;  therefore,  the  gyro  can  be  made  to  tumble 
about  two  axes.  In  Figure  1,  SRA  and  OA  indicate  direction  of  spin 
reference  axis  and  output  axis,  respectively. 

The  test  is  performed  with  the  table  axis  positioned  parallel  to 
the  earth  rotational  axis.  The  first  table  axis  is  elevated  to  the 
local  latitude  angle  A from  the  horizontal  plane.  Under  this  condition, 
IA  is  perpendicular  to  polar  axis  and  the  gyro  will  not  sense  the  earth 
rata. 

Let  the  first  table  be  rotated  at  an  angular  rate  to,  and  the  second 
table  be  rotated  about  the  IA  axis  at  an  angular  rate  to  . The  gyro 


senses  o>  but  not  ox  . Let  ox,  be  the  gyro  output.  Because  to,  is  known, 

0 <p  K 0 

the  difference  between  toD  and  to  , which  is  the  total  drift,  can  be 
calculated.  In  equation  form. 


D = (Or  - (O0 


(5) 


To  avoid  the  earth  rate  effect  which  may  be  caused  by  imperfect  align- 
ment of  the  table  axis  to  the  polar  axis,  it  is  desirable  to  turn  the 
table  at  a rate  approximately  ten  to  twenty  times  faster  than  the  earth 
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rate.  To  have  a good  frequency  separation  property  for  the  drift 
signal,  which  will  be  discussed  later,  the  rotation  rate  of  the  second 
axis  should  be  approximately  five  times  that  of  the  table  axis.  A 
reasonable  data  taking  interval  is  every  fet0<feLO  the  second  axis  ^ 

rotation.  More  data  points  allow  a better  statistical  data  reduction. 
The  required  data  record  covers  a full  rotation  of  the  first  axis,  or 
one  of  its  integer  multiples. 

The  remaining  problem  is  the  determination  of  nine  drift  parameters 
from  the  measured  values  of  D by  means  of  computer  data  processing. 


IV.  DATA  PROCESSING  ALGORITHM 

A.  The  Measurement  Equation 

Equation  (4)  shows  that  eight  of  the  nine  drift  terms  are 
dependent  on  either  components  of  acceleration  or  on  their  produc.s. 
Because  the  test  stand  is  stationary,  the  acceleration  components  experi- 
enced by  the  gyro  are  due  to  earth  gravitation.  The  value  of  each  accel- 
eration component,  in  turn,  depends  on  the  angular  positions  of  the  first 
and  second  table  axes.  These  angular  positions  are  represented  by 
$ and  0 as  shown  in  Figure  1 where  their  references  of  measurement  are 
also  indicated. 

The  acceleration  components  experienced  by  the  gyro  can  be  obtained 
with  the  help  of  Figure  2.  The  accelerations  along  the  spin  reference 
axis,  input  axis,  and  output  axis  are  given,  respectively,  by 


a^  = g cos  A cos  0 cos  6 + g sin  A sin  0 


a^  = -g  cos  A sin  <f> 


aQ  - g cos  A cos  <f>  sin  0 - g sin  A cos  B 


(6) 

(7) 

(8) 


where  g is  the  earth  gravitation.  Let 


0 = m0t 


<*>  = 

<P 


(9) 

(10) 
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Then  Equations  (6),  (7),  and  (8)  become 


a = g cos  X cos  ioQt  cos  u.t  + g sin  X sin  w.t 
s 09  o 


a.  = -g  cos  X sin  <o,t 

x 9 


a = g cos  X cos  to,t  sin  to.t  - g sin  X cos  to.t 
o 90  o 


Substituting  the  three  preceding  equations  into  Equation  (4)  gives 
the  total  drift  as  a function  of  time  t;  or,  equivalently  as  a function 
of  the  angular  positions  of  the  table  and  the  gyro  house.  Thie  explicit 
.form  of  the  total  drift  is 


D = D + K {g  cos  X cos  to.t  cos  co,t  + g sin  X sin  to,t} 
os  6 6 c 


- {g  cos  X sin  <o^t) 


+ K {g  cos  X cos  w.t  sin  u>  t - g sin  X cos  mfit} 


+ K {g  cos  X cos  to.t  cos  to.t  + g sin  X sin  io_t}' 
ss  09  v 


+ K..  {g  cos  X sin  m.t}' 
ix  9 


K . {g  cos  X cos  u)„t  cos  to  t + g sin  X sir  u.t}  x {g  cos  X sin  u t} 

SI  0 9 5 9 

K.  {g  cos  X sin  wxt}  x {g  cos  X cos  u t sin  tc.t  - g sin  X cos  a.t} 
10  9 90 


K {g  sin  X cos  wnt  - g cos  X cos  to  t sin  to  t} 
os  0 9 o 


x {g  cos  X cos  to  t cos  to  t + g sin  X sin  m.t} 


Equation  (14)  contains  terms  which  are  products  of  two  frequencies, 
indicating  the  existance  of  modulation  process.  As  a result,  extra 
frequency  components  are  generated.  With  the  aid  of  trigonometric 
identities.  Equation  (14)  can  be  rearranged  as  a sum  of  several  groups; 
each  is  associated  with  one  frequency  as  follows; 
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Equation  (15)  is  the  measurement  equation  sought.  The  equation  relaces 
the  total  drift  D to  drift  parameters  which  are  Dq  and  K's.  It  is  noted 

that  the  total  drift  contains  thirteen  different  frequency  components. 

The  right-hand  side  of  the  equation  is  arranged  into  thirteen  terms;  each 
term  gives  one  frequency  component. 

B.  Use  of  Fourier  Analysis 

If  values  of  <o,  and  are  property  chosen,  the  thirteen 

<p  v 

frequencies  in  the  drift  signal  can  be  distinctly  separated  from  each 
other.  Then  the  amplitudes  and  phases  of  respective  frequency  components 
can  be  determined  from  the  measured  drift  signal  by  performing  a Fourier 
analysis  of  the  signal.  This  can  be  effected  by  using  either  a Fast 
Fourier  transformation  (FFT)  or  by  Fourier  coefficient  determination  [8]. 
In  terms  of  amplitudes  and  phases.  Equation  (15)  has  the  form 


D = A.  + A0  sin  (w  t + B0)  + A_ 
12  <p  2 o 


+ A^,  sin  [(u>0  - 2w^)  t + B^]  + A5  sin  [ (wQ  - w^)t  + 


+ A6  sin  (w0t  + B0)  + A?  sin  [ (w0  + o^)t  + 87] 


0 r 


+ Ag  sin  [(u)0  + 2(0^)  t + 3g]  + A$  sin  [2(w0  - w^t  + 8g] 


+ A1q  sin  [ (2m0  - w^t  + &1Ql  + sin  (2u>0t  + B^) 


+ A^^  sin  t(2w0  + u^)  + B12J  + A^  sin  [2(w0  + w^)t  + 3^]  (17) 


where 


A = D + K 

1 O S£ 


A2  = KiD4 


+ D2/+KiiDl 


A3  = Kss  T - KiiDl 


EIS^ 


e6  " t“'1  { 


M 

K / 

o • 


(4th  q^) 


(36) 


-1  Ks°4  + Kai°3 

e7  tan  K D + K D,  (1St  ql} 
o 4 io  3 


(37) 


h" 

■ tan"1 ( 

io 

K 

si  ' 

(2nd  qj) 

(38) 

V 

. tan’1  ( 

K \ 
ss  j 

' os ' 

(1st  q2) 

(39) 

eio 

„ "I 

= tan 

ol  03 

* 

(4th  qj 

(40) 

ell 

= tan  ^ 

ft) 

(1st  qx) 

(41) 

ei2 

= tan  ^ 

(ft) 

ss  ' 

(4th  qx) 

(4-) 

813 

= tan  ^ 

ft) 

os 

(1st  q1) 

(43) 

The  parentheses  following  each  angle  indicate  the  quadrant  of  the  angle. 

Equations  (18)  to  (43)  consist  of  a total  of  twenty-six  equations. 
They  are  used  to  solve  for  the  nine  drift  coefficients  after  and  B^, 

i = 1,  ...,  13,  are  determined  from  the  Fourier  analysis  of  the  drift 
signal.  There  are  more  equations  than  unknowns,  indicating  redundant 
information.  Maximum  use  of  the  information  will  minimize  the  effect 
of  instrument  noise  on  the  determination  of  drift  coefficients  from  the 
Fourier  coefficients.  The  method  chosen  for  reducing  the  redundant 
information  to  an  unique  set  of  drift  coefficients  is  the  well-known 
least  square  regression  [9,  10]. 

C.  Use  of  Least  Square  Regression 

Squaring  both  sides  of  Equations  (26)  to  (30),  combining 
the  result  with  Equations  (39)  to  (43),  and  arranging  these  ten  equations 
in  a matrix  equation  form,  gives 
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± 


u M 


Vector  u_  and  matrix  M are  defined  as  shown  in  Equation  (44).  A well- 
known  least  square  regression  formula  gives  [9] 


= (MTM)-1  MTu  . (45) 


By  taking  the  square  root  of  Equation  (45)  and  retaining  the  positive 

values,  K and  K are  determined. 

’ ss  os 

Similarly,  squaring,  summing,  and  rearranging,  Equations  (21)  to 
(25)  and  (34)  to  (38)  can  give 


where 

a = tan  8^  tan  6^  . (47) 

Vector  v and  matrix  N are  defined  as  shown  in  Equation  (46) ; the  least 
square  solution  for  the  equation  is 


Taking  the  positive  values  of  the  square  root  of  Equation  (48)  gives 
K v K . , and  K . 


V 


Finally,  using  Equations  (20),  (19),  and  (18)  in  this  order,  the 
following  formulas  are  obtained: 


" Kii°l 


(49) 


(50) 

(51) 


Thus,  all  nine  drift  parameters  are  determined. 

D.  Discussion 

A striking  difference  between  the  new  and  the  old  two- 
axis  tumble  tests  is  that  the  former  does  not  require  the  knowledge  of 
the  time  instant  at  each  data  point  (or  equivalently,  the  attitude  of 
the  second  axis  at  each  data  point).  This  allows  the  use  of  simpler 
test  equipment  and  demands  less  operator  effort.  The  conceptual  differ- 
ence is  that  in  the  new  method  the  drift  parameter  identification  is 
performed  using  the  frequency  domain  information  while  in  the  latter 
method  identification  is  done  using  the  time  domain  information.  It 
should  be  mentioned  that  the  new  method  reduces  the  operator's  effort 
and  lessens  the  test  equipment  requirement  at  the  expense  of  a slightly 
more  complex  data  processing  algorithm.  Because  both  the  new  and  the 
old  method  require  a computer  for  data  processing  and  because  the  computer 
time  required  for  either  method  is  insignificant,  the  trade-off  favors 
the  new  method. 


On  the  other  hand,  both  the  new  and  old  two-axis  tumble  tests  have 
the  following  advantages  as  compared  to  the  other  three  methods  listed 
in  the  Introduction.  First,  there  is  no  need  of  repeated  gyro  reposi- 
tioning which  is  time-consuming.  Second,  because  both  the  first  and 
second  tables  are  turning  at  constant  speed,  there  is  no  motion-induced 
transient  effect  in  the  measured  drift  signal.  Third,  the  measured 
drift  signal  provides  sufficient  information  for  the  determination  of 
all  nine,  rather  than  partial,  drift  parameters.  Fourth,  the  required 
test  time  is  shorter. 


Finally,  to  be  truly  objective,  it  should  be  mentioned  chat  the  old 
two-axis  tumble  test  has  one  advantage  over  the  new  method.  When  using 
the  old  method,  the  required  data  record  can  be  of  any  length  as  long  as 
it  contains  at  least  nine  data  points.  In  other  words,  the  test  operator 
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does  not  have  to  wait  for  the  table  to  turn  a full  revolution  although 
the  reduced  data  record  may  not  provide  sufficient  redundancy  for  a 
good  statistical  data  processing.  The  algorithm  for  the  old  method  is 
included  as  the  Appendix  for  the  reader's  reference. 


V.  SIMULATION  RESULT 

Simulated  testing  was  conducted  to  demonstrate  the  proposed 
method.  The  simulation  wao  performed  with  the  first  table  turning 
about  its  axis  at  twenty  times  the  earth  rate  and  the  second  table 
turning  about  its  axis  at  100  times  the  earth  rate.  The  gyro  drift 
was  simulated.  For  a full  rotation  of  the  table,  864  data  points  were 
generated. 

With  the  chosen  rotation  rates  for  the  two  tables,  the  generated 
frequency  components  in  the  drift  signal  are  evenly  separated  by  a 
frequency  separation  of  20  Hz.  Table  1 lists  all  the  frequency  components. 
The  amplitudes  and  tangent  of  the  phases  of  the  frequency  components  as 
generated  by  Fourier  analysis  are  shown  in  Table  2.  All  nine  drift 
coefficients  resulting  from  the  least  square  regression  are  shown  in 
Table  3.  The  effectiveness  of  the  identification  technique  and  the 
correctness  of  the  data  reduction  algorithm  are  demonstrated.  The 
discrepancies  are  evidently  due  to  computation  errors,  which  can  be 
reduced  by  using  the  double  precision  computation.  However,  the  identi- 
fication accuracy  as  shown  in  Table  3 is  sufficient  for  most  applications. 


VI.  CONCLUSION 


A frequency  domain  technique  for  the  identification  of  gyro 
drift  coefficients  has  been  developed.  The  technique  involves  a 
two-axis  tumble  test  and  an  unique  data  reduction.  Compared  to  other 
techniques,  this  technique  has  the  merit  of  providing  more  complete 
data  for  describing  gyro  drift  characteristics.  In  addition,  the 
required  procedure  for  testing  and  data  reduction  can  be  conveniently 
implemented  into  an  automated  gyro  drift  coefficient  identification 
system. 

A simulation  was  performed.  The  result  demonstrated  the  effective- 
ness of  the  gyro  drift  coefficient  determination  technique  and  the 
correctness  of  the  data  reduction  algorithm. 
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TABLE  1.  FREQUENCY  COMPONENTS  IN  IRIFT  SIGNAL 


Frequency 

Component 

Frequency 

Value 

(Hz) 

d-c 

0 

“♦ 

20 

% 

40 

W8  ~ % 

60 

-&■ 

3 

1 

<D 

3 

80 

we 

100 

“e  + «* 

120 

“e  + \ 

140 

2(aj0  ’ V 

160 

2u9  ~ % 

180 

2ue 

200 

2we  + "$ 

220 

2(w„  +U,) 

240 

TABLE  2.  RESULT  OF  FOURIER  ANALYSIS 


Frequency 

Real  Part 

Imaginary 

Part 

de 

0.82279050 



20 

-0.00000126 

-0.22565949 

40 

0.00138369 

-0.00000834 

60 

-0.00023610 

0.00023035 

80 

0.08428782 

0.01377910 

100 

-0.01948988 

0.11692071 

120 

0.08493614 

0.01442608 

140 

0.00023283 

-0.00023330 

160 

0.00023244 

0.00011725 

180 

-0.00032030 

0.00064254 

200 

0.00002115 

0.00001165 

220 

-0. 00032045 

0.00064246 

240 

0.00023199 

0.00011721 

20 


^'35fe^g£%33^ 


TABLE  3.  IDENTIFIED  DRIFT  COEFFICIENTS 


Drift 

Coefficient 

Value 

Discrepancy 

(%) 

Identified 

.... ...  ... 

Actual 

Dq  (deg/hr) 

5. 99982C 

j 

6.0 

0.003 

Kg  (deg/hr/g) 

2.999528 

3.0 

0.016 

Ki  (deg/hr/g) 

4.000009 

4.0 

0.0002 

Kq  (deg/hr/g) 

0.5000007 

0.5 

0.0001 

Kgg  (deg/hr/g2) 

0.040064 

0.04 

0.16 

KijL  (deg/hr/g2) 

-0.039625 

-0.04 

0.94 

Kgi  (deg/hr/g2) 

0.020106 

0.02 

0.53 

Kio  (deg/hr/g2) 

0.0200658 

0.02 

0.33 

K (deg/hr/g2) 

OS 

0.0202097 

0.02 

1.05 
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Appendix.  ALGORITHM  FOR  OLD  TWO-AXIS  TUMBLE  TEST 

Let  a and  8 be  the  incremental  angular  displacements  of  the  table 
and  gyro  house,  respectively,  between  two  consecutive  data  taking  points. 
Let  <£  = 8 = 0 initially.  Then  at  the  k-th  data  taking  point,  <£  = ke 
and  6 = k8.  Recalling  the  definition  of  and  in  Equation  (16),  then 

Equations  (6),  (7),  and  (8)  can  be  expressed  as 


a (k)  = D,  cos  ku  cos  kS  + D_  sin  k8 
s 4 5 


a_^(k)  = -D^  sin  ka 


a (k)  = D,  cos  ka  sin  k8  - Dc  cos  k8 
o 4 5 


Define 


bss(k)  = ag(k) 


bii.(k)  = ai^k) 


b . (k)  = a (k)  a.  (k) 
si  s 1 


b.  (k)  = a.(k)  a (k) 
10  i o 


bos(k)  = ao(k)  as(k) 


Substituting  Equations  (55)  through  (59)  into  Equation  (4)  gives 


D(k)  = D + K a (k)  + K.  a. (k)  + K a (k) 
' o s s 11  o o' 


+ K b (k)  + K. . b. . (k) 
ss  ss  11  11 


+ K . t .(k)  + K.  b.  (k)  + K b (k) 

SI  SI  lO  10  OS  OS 


Using  vector  notation. 


D(k)  = CT(k)  x 


(61) 


where 


C(k)  = 


as(k) 

ai(k) 

ao(k; 

bss<k> 

b.100 

b .(k) 
sxv 

b.  (k) 
xo 

b (k) 
os 


(62) 


x = 


K 


K. 

i 

K 


K 

ss 

K. . 

xx 

K . 
sx 

K. 

xo 


K 

os 


(63) 


In  Equation  (61),  D(k)  is  the  measured  drift  at  the  k-th  time 
instant  and  x is  the  parameter  vector  to  be  determined.  For  a total  of 
n measurements,  there  are  n equations  in  the  form  of  Equation  (61).  They 
can  be  arranged  in  a column  as  follows: 
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• T 

C (1) 


D(k)  “ CA(k)  x 


D(n)J  U (n). 


Define  vector  _z  and  matrix  Q as  shown  in  Equation  (64).  In  general, 
n > 9,  so  Q is  not  square  and  therefore  cannot  be  inverted.  An  unique 
solution  for  x can  be  obtained  using  the  pseudo  inverse 


x = (QTQ)_1QTz 


which  is  exactly  the  same  p. o the  least  square  regression  formula  used  in 
Equations  (45)  and  (48) . 

Equation  (65)  is  a batch  processing  algorithm.  It  involves  multi- 
plication of  n x 9 matrices  and  inversion  of  a 9 * 9 matrix.  The 
algorithm  can  be  made  sequential  if  so  desired. 
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